Creating maps in R - basics
Heidi Rautiainen
2024-12-20, Updated: 2026-01-21
Set up directory, install packages and load packages
Exercise 1: Create study area map using “ggplot”
All datasets are found in the “data”-folder.
Here, we use the package ‘rnaturalearth’ to download countries in a ESRI shapefile format. These are available for download here: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/
Coordinate reference systems (CRS)
CRS provide a standardized way of describing locations to describe
geographic data. The CRS that is chosen depends on when the data was
collected, the purpose of the data, etc. It is necessary to transform
vector and raster data to a common CRS when data with different CRS are
combined.
Setting CRS to SWEREF99 (Sweden). CRS can be referenced by its:
EPSG code, CRS(“+init=epsg:3006”) (see http://www.epsg-registry.org/ and http://spatialreference.org/) or by
proj4spring; “+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”
There are two general options:
unprojected (a.k.a. Geographic): Latitude/Longitude for referencing location on the ellipsoid Earth (wgs84), and
projected: Easting/Northing for referencing location on 2D representations of Earth (the creation of maps) e.g., SWEREF99.
More reading:
https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf
Load vector data
How to download using “rnaturalearth” and save the shp locally
I will open the saved .shp for the exercise:
## load country boarders -----------
#' downladed and saved from rnaturalearth (above)
swe_no_sweref <- input_vector %>%
# value = T: return vector containing the matching elements
grep(pattern = "swe_no", value = T) %>%
# read ESRI Shapefile object
sf::st_read() %>%
#convert geometry object into an sf object
st_as_sf() ## Reading layer `swe_no' from data source
## `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/country/swe_no.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -362704.5 ymin: -6097519 xmax: 1114145 ymax: 8972413
## Projected CRS: SWEREF99 TM# check coordinate system
st_crs(swe_no_sweref)$proj4string## [1] "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"# Create a bounding box as an sf object for layers ----
b_box <- c(xmin = 350432, ymin = 7179284,
xmax = 857780, ymax =7444230)
# Now intersect the full world's coastline with the bounding box
suppressWarnings({
sweden.c <- st_crop(swe_no_sweref, b_box)
plot(sweden.c)
}) 
class(sweden.c)## [1] "sf" "data.frame"### water data ----
input_vector## [1] "input_data/vector/animal_data/data_cleaned_PV.rds"
## [2] "input_data/vector/animal_data/data_ready_rsf.rds"
## [3] "input_data/vector/animal_data/reindeer_data.rds"
## [4] "input_data/vector/country/swe_no.shp"
## [5] "input_data/vector/country/swe.shp"
## [6] "input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp"
## [7] "input_data/vector/ne_10m_lakes/ne_10m_lakes.shp"# import lakes
lakes <- input_vector %>%
grep(pattern = "lakes", value = T) %>%
sf::st_read() %>%
st_as_sf() ## Reading layer `ne_10m_lakes' from data source
## `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/ne_10m_lakes/ne_10m_lakes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1355 features and 41 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -165.9656 ymin: -50.66967 xmax: 177.1544 ymax: 81.95521
## Geodetic CRS: WGS 84class(lakes)## [1] "sf" "data.frame"st_crs(lakes)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]# Transform or convert coordinates of simple feature (sf)
lakes <- st_transform(lakes, crs = 3006)
suppressWarnings({
# crop an sf object to a specific rectangle (here based on extent of "sweden.c"-shp)
# sf_use_s2(FALSE) # use if errors occur below
lakes <- st_crop(lakes, sweden.c)
class(lakes)
plot(lakes)
}) 
Create points of interest (cities)
Create xy points with labels. Here creating sf-object and setting the crs.
# Create points of interest (coordinates from QGIS or google maps for demonstration)
places <- data.frame(ID= c("Arjeplog", "Arvidsjaur", "Sorsele"),
y = c(7328407.945, 7280856.153, 7270326.684),
x = c(630040.101, 693262.414, 617477.638)) %>%
st_as_sf(coords = c("x", "y"), # create sf object
crs=st_crs(sweden.c)) # set same crs as for sweden.c
class(places)## [1] "sf" "data.frame" st_crs(places) # check projection ## Coordinate Reference System:
## User input: SWEREF99 TM
## wkt:
## PROJCRS["SWEREF99 TM",
## BASEGEOGCRS["SWEREF99",
## DATUM["SWEREF99",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4619]],
## CONVERSION["SWEREF99 TM",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Topographic mapping (medium and small scale)."],
## AREA["Sweden - onshore and offshore."],
## BBOX[54.96,10.03,69.07,24.17]],
## ID["EPSG",3006]]Make plot using ggplot
library(ggspatial)
plot(sweden.c)## Warning: plotting the first 9 out of 168 attributes; use max.plot = 168 to plot
## all
plot(lakes)## Warning: plotting the first 9 out of 41 attributes; use max.plot = 41 to plot
## all
plot(places)
ggplot() +
# plot sweden and lakes
geom_sf(aes(geometry = sweden.c$geometry)) + # note: geom_sf() for vector data
geom_sf(data = sweden.c, # plot sweden border vector
aes(geometry = sweden.c$geometry),
fill= "antiquewhite")+
geom_sf(data = lakes, # plot lakes vector
fill="lightblue") +
# layer_spatial(data=elevation) +
# add cities
geom_sf(data = places) +
# add city names
geom_sf_text(data = places, aes(label = ID),
vjust=-0.5)+
# Add text on map
annotate(geom = "text",
x = 484609, y = 7404230,
label = "Norway",
fontface = "italic",
color = "darkgray",
size = 6) +
annotate(geom = "text",
x = 784609, y = 7404230,
label = "Sweden",
fontface = "italic",
color = "darkgray",
size = 6) +
# Add north arrow
annotation_north_arrow(location = "bl",
which_north = "true",
height = unit(1, "cm"),
width = unit(1, "cm"),
pad_x = unit(0.5, "cm"), # horizontal align
pad_y = unit(8, "cm"), # vertical align
style = north_arrow_fancy_orienteering) +
#add annotation scale
annotation_scale(location = "bl",
width_hint = 0.5) +
# keep new projection for ggplot
coord_sf(crs = st_crs(3006),
datum = sf::st_crs(3006),
expand = FALSE) +
# add grid lines
theme(panel.grid.major = element_line(color = gray(.5),
linetype = "dashed",
linewidth = 0.1),
panel.background = element_rect(fill = "lightblue")) 
### examplefor removing or adding titles
# theme(axis.title.x = element_blank(), # remove
# axis.title.y = element_blank()) +
# xlab("Longitude") + ylab("Latitude") # add
ggsave("output_data/Fig1_study_area.png",
width = 8,
height = 6,
dpi = 400)Optional: Create a polygon and add to map
Example 2 interactive maps using mapview
# some gps data from reindeer
reindeer_data <- readRDS("input_data/vector/animal_data/reindeer_data.rds") %>%
dplyr::select(Collar_ID, id, year, x_, y_, t_, group)
str(reindeer_data)## 'data.frame': 5951 obs. of 7 variables:
## $ Collar_ID: chr "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
## $ id : chr "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
## $ year : int 2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
## $ x_ : num 733706 733653 733836 733828 733683 ...
## $ y_ : num 7314782 7314737 7314807 7314851 7314738 ...
## $ t_ : chr "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
## $ group : chr "Fed" "Fed" "Fed" "Fed" ...### basic usage ===========================================================
### most basic call
mapview()### with some data
# mapview(reindeer_data)
str(reindeer_data)## 'data.frame': 5951 obs. of 7 variables:
## $ Collar_ID: chr "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
## $ id : chr "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
## $ year : int 2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
## $ x_ : num 733706 733653 733836 733828 733683 ...
## $ y_ : num 7314782 7314737 7314807 7314851 7314738 ...
## $ t_ : chr "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
## $ group : chr "Fed" "Fed" "Fed" "Fed" ...# convert to sf
reindeer_data_sf <- reindeer_data %>%
st_as_sf(coords = c("x_", "y_"), #data frame to sf
crs=3006) # assign to known crs
plot(reindeer_data_sf)
mapview(reindeer_data_sf) # mapview structure
m <- mapview(reindeer_data_sf)
mstr(m, 3) # two slots = object (points) and map (geometries)## Formal class 'mapview' [package "mapview"] with 2 slots
## ..@ object:List of 1
## .. ..$ :sfc_POINT of length 5951; first list element: 'XY' num [1:2] 20.1 65.9
## ..@ map :List of 8
## .. ..$ x :List of 4
## .. ..$ width : NULL
## .. ..$ height : NULL
## .. ..$ sizingPolicy :List of 7
## .. ..$ dependencies :List of 5
## .. ..$ elementId : NULL
## .. ..$ preRenderHook:function (widget)
## .. ..$ jsHooks :List of 1
## .. ..- attr(*, "class")= chr [1:2] "leaflet" "htmlwidget"
## .. ..- attr(*, "package")= chr "leaflet"m@object## [[1]]
## Geometry set for 5951 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 19.75501 ymin: 65.85719 xmax: 20.24619 ymax: 65.92977
## Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
## First 5 geometries:## POINT (20.12738 65.86829)## POINT (20.12616 65.86794)## POINT (20.13028 65.86842)## POINT (20.13018 65.86882)## POINT (20.12679 65.86792)class(m@map)## [1] "leaflet" "htmlwidget"# mapview
mapview(sweden.c) # polygons# character to factor needed for grouping colors
reindeer_data$group <- as.factor(reindeer_data$group)
unique(reindeer_data$Collar_ID)## [1] "FIT_421952" "FIT_421943" "FIT_421941" "FIT_421940" "FIT_421949"
## [6] "FIT_421910" "FIT_421939" "FIT_421946" "FIT_421906" "FIT_421956"
## [11] "FIT_421944" "FIT_421937" "FIT_421948" "FIT_421908" "FIT_421935"
## [16] "FIT_421942" "FIT_421954" "FIT_421907" "FIT_421932" "FIT_421931"
## [21] "FIT_421934" "FIT_421945" "FIT_421938" "FIT_421947" "FIT_421936"
## [26] "FIT_421953" "FIT_421909" "FIT_421950" "FIT_421951" "FIT_421933"library(mapview)
mapid <- reindeer_data %>%
filter(Collar_ID =="FIT_421952" | Collar_ID =="FIT_421943"| Collar_ID =="FIT_421941") %>%
st_as_sf(coords = c("x_", "y_"), crs=3006)
str(mapid)## Classes 'sf' and 'data.frame': 628 obs. of 6 variables:
## $ Collar_ID: chr "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
## $ id : chr "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
## $ year : int 2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
## $ t_ : chr "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
## $ group : Factor w/ 2 levels "Control","Fed": 2 2 2 2 2 2 2 2 2 2 ...
## $ geometry :sfc_POINT of length 628; first list element: 'XY' num 733706 7314782
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
## ..- attr(*, "names")= chr [1:5] "Collar_ID" "id" "year" "t_" ...### Add home range
AOI <- grep("modified_new.shp", input_vector, value = T) %>%
sf::read_sf() %>%
sf::st_transform(crs = 3006)
plot(AOI)## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
pal <- mapviewPalette("mapviewVectorColors") # optional: mapviewTopoColors, mapviewRasterColors, mapviewVectorColors, mapviewSpectralColors
mapid %>%
mapview(zcol="Collar_ID",
cex = "Collar_ID",
col.regions = pal(3))+
mapview(
AOI,
alpha.regions = 0, # no fill
color = "black", # border color
lwd = 2 # border width
) # or simply add + AOI ### view extents
viewExtent(mapid)mapview(st_bbox(mapid))viewExtent(mapid) + mapview(mapid, zcol = "Collar_ID")#### show tracks
mapid2 <- reindeer_data %>%
filter(Collar_ID %in% c("FIT_421952","FIT_421943","FIT_421941")) %>%
mutate(t_ = as.POSIXct(t_, format = "%Y-%m-%d %H:%M:%S", tz = "UTC")) %>%
arrange(Collar_ID, t_) %>%
st_as_sf(coords = c("x_", "y_"), crs = 3006)
str(mapid2)## Classes 'sf' and 'data.frame': 628 obs. of 6 variables:
## $ Collar_ID: chr "FIT_421941" "FIT_421941" "FIT_421941" "FIT_421941" ...
## $ id : chr "rt_PV_20_024" "rt_PV_20_024" "rt_PV_20_024" "rt_PV_20_024" ...
## $ year : int 2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
## $ t_ : POSIXct, format: "2021-01-29 20:00:44" "2021-01-30 00:01:20" ...
## $ group : Factor w/ 2 levels "Control","Fed": 2 2 2 2 2 2 2 2 2 2 ...
## $ geometry :sfc_POINT of length 628; first list element: 'XY' num 733461 7315029
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
## ..- attr(*, "names")= chr [1:5] "Collar_ID" "id" "year" "t_" ...tracks <- mapid2 %>%
group_by(Collar_ID) %>%
summarise(do_union = FALSE) %>% # keeps point order within groups
st_cast("LINESTRING")
mapview(tracks,
zcol = "Collar_ID",
lwd = 3,
layer.name = "Tracks") +
mapview(mapid, zcol = "Collar_ID", cex = 3, layer.name = "Fixes") ### basemaps
library(leaflet)
# To specify which background map you would like to use
mapview(gadmCHE, map.types = c("Stamen.Toner", "NASAGIBS.ViirsEarthAtNight2012"))# see https://leaflet-extras.github.io/leaflet-providers/preview/ for more options -- to see which background maps that leaflet provides! Exercise 3 (optional): Maps in R
Pre-define crs using ESPG-code or proj4string by:
ESPG: CRS(“+init=epsg:3006”)
proj4string: “+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs”
# projection to use as proj4string
crs_to_use <- "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs" Load background data
Now let’s start by loading the maps using the “raster”-package. Another option: “terra”-package using rast()-function)
## [1] "input_data/ignore/PV_2021_cc.tif"
## [2] "input_data/ignore/PV_2021_houses.tif"
## [3] "input_data/ignore/PV_2021_lichen.tif"
## [4] "input_data/ignore/PV_2021_nmd.tif"
## [5] "input_data/ignore/PV_2021_roads.tif"
## [6] "input_data/ignore/PV_2021_TRI.tif"
## [7] "input_data/raster/PV_2021_dem.tif"
## [8] "input_data/raster/PV_2021_slope.tif"## [1] "Clear_cuts_2021" "prox10_houses_stakke"
## [3] "lichenmap_stakke_original" "nmd10_stakke_ungeneralized"
## [5] "prox10_stakke_resample" "ruggedness10_QGIS"
## [7] "PV_2021_dem" "slope10m_QGIS"## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"
## Reading layer `spatial_files_availability_PV_winter_2021_shorter_modified_new' from data source `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 716285.4 ymin: 7313268 xmax: 739293.8 ymax: 7321582
## Projected CRS: SWEREF99 TM## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"
Add animal data (GPS points)
Load prepared reindeer data
## 'data.frame': 5951 obs. of 7 variables:
## $ Collar_ID: chr "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
## $ id : chr "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
## $ year : int 2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
## $ x_ : num 733706 733653 733836 733828 733683 ...
## $ y_ : num 7314782 7314737 7314807 7314851 7314738 ...
## $ t_ : chr "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
## $ group : chr "Fed" "Fed" "Fed" "Fed" ...## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
Explore and select colors
## [1] "#00A600" "#01A600" "#03A700" "#04A700" "#05A800" "#07A800"
Plot and save all together
## Warning in plot.sf(HR_area, col = NA, border = 1, add = T): ignoring all but
## the first attribute
Some optional adjustments: hillshade instead of elevation
Prepare hillshade background map
Calculate slope and aspect from elevation (dem)
Plot and save all together
## Warning in plot.sf(HR_area, col = NA, border = 1, add = T): ignoring all but
## the first attribute